#include "trans_coeff.h"
#include "nvtx3/nvToolsExt.h"
#include <cmath>

int TransCoeff::phase(int n) { return n % 2 ? -1 : 1; }
int TransCoeff::Tri(int &j2a, int &j2b, int &j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}
int TransCoeff::Tri(int j2a, int j2b, int j2c) {
  // I+J>=K, I+K>=J, J+K>=I,
  // I/2+J/2+K/2 = INTEGER.
  // TRI=1, WHEN TRIADIC CONDITION IS FULFILLED, TRI=-1 OTHERWISE
  int L2 = j2a + j2b + j2c;
  if (int(L2) != (L2 / 2) * 2)
    return -1;
  L2 = L2 / 2;
  if (j2a * j2b * j2c < 0)
    return -1;
  if ((L2 - j2a) * (L2 - j2b) * (L2 - j2c) < 0)
    return -1;
  return 1;
}

double TransCoeff::cal(State_Jacobi &jaco, State_TP &triP) {
  double result = 0.0;
  if (jaco.T_2 != triP.T_2 || jaco.J2_NL != triP.J_2)
    return 0.0;

  if (jaco.E_NL != triP.E)
    return 0.0;

  if (jaco.T12_2 != triP.Tab_2)
    return 0.0;
  // cout<<"good"<<endl;

  int a_x = triP.a;
  int b_x = triP.b;
  int c_x = triP.c;
  int l_a = SP.state[a_x].l;
  int l_b = SP.state[b_x].l;
  int l_c = SP.state[c_x].l;
  int E_a = SP.state[a_x].E;
  int E_b = SP.state[b_x].E;
  int E_c = SP.state[c_x].E;
  int Eab = SP.state[a_x].E + SP.state[b_x].E;
  int E12 = 2 * jaco.n12 + jaco.l12;
  int E3 = 2 * jaco.n3 + jaco.l3;
  int E0 = 2 * jaco.N_0 + jaco.L_0;

  // cout<<"a :"<<a_x<<endl;
  // cout<<"b :"<<b_x<<endl;
  // cout<<"c :"<<c_x<<endl;

  double d1 = 1.0;
  double d2 = 2.0;

  int L1c_min = std::abs(jaco.L_0 - jaco.l3);
  int L1c_max = jaco.L_0 + jaco.l3;

  int L123_min = std::abs(jaco.l12 - jaco.l3);
  int L123_max = jaco.l12 + jaco.l3;

  // part1
  double part1 = (SP.state[a_x].j_2 + 1) * (SP.state[b_x].j_2 + 1) *
                 (SP.state[c_x].j_2 + 1) * (triP.Jab_2 + 1);
  part1 = part1 * (jaco.J_2 + 1) * (jaco.J12_2 + 1) * (jaco.j3_2 + 1);
  part1 = std::sqrt(part1);
  int phase_x = SP.state[c_x].l + jaco.l12 + jaco.L_0;
  part1 = part1 * phase(phase_x);
  // cout<<"part1 = "<<part1<<endl;
  // part2 part3 ...
  int Lab_min = std::abs(l_a - l_b);
  int Lab_max = l_a + l_b;

  // cout<<"\t"<<L1c_min<<"\t"<<L1c_max;
  // cout<<"\t"<<L123_min<<"\t"<<L123_max;
  // cout<<"\t"<<Lab_min<<"\t"<<Lab_max<<endl;

  double part2 = 0.0;
  for (int Sab_2 = 0; Sab_2 <= 2; Sab_2 = Sab_2 + 2) {
    int S2_min = std::abs(Sab_2 - 1);
    int S2_max = Sab_2 + 1;
    for (int S_2 = S2_min; S_2 <= S2_max; S_2 = S_2 + 2) {
      for (int Lab = Lab_min; Lab <= Lab_max; Lab++) {
        int L_min = std::abs(Lab - l_c);
        int L_max = Lab + l_c;
        // cout<<"L : "<<L_min<<"\t"<<L_max<<endl;
        for (int L = L_min; L <= L_max; L++) {
          double part2_t = 0.0;
          part2_t = ninJ.cal(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, Sab_2,
                             SP.state[a_x].j_2, SP.state[b_x].j_2, triP.Jab_2);
          // cout<<"A = "<<part2_t<<endl;
          if (part2_t == 0.0)
            continue;
          part2_t = part2_t * ninJ.cal(2 * Lab, 2 * l_c, 2 * L, Sab_2, 1, S_2,
                                       triP.Jab_2, SP.state[c_x].j_2, triP.J_2);
          // cout<<Lab<<"\t"<<l_c<<"\t"<<L<<"\t"<<Sab_2<<"\t"<<1<<"\t"<<S_2<<"\t";
          // cout<<triP.Jab_2<<"\t"<<SP.state[c_x].j_2<<"\t"<<triP.J_2<<endl;
          // cout<<"A*B = "<<part2_t<<endl;
          if (part2_t == 0.0)
            continue;
          part2_t = part2_t * (2 * Lab + 1) * (2 * L + 1) * (S_2 + 1) *
                    std::sqrt(Sab_2 + 1);
          phase_x = Lab + L + (S_2 + jaco.J_2) / 2;
          part2_t = part2_t * phase(phase_x);

          if (part2_t != 0.0) {
            // cout<<"check point 1"<<endl;
          }
          int L1_min = std::abs(Lab - jaco.l12);
          int L1_max = Lab + jaco.l12;
          // part3
          double part3 = 0.0;
          for (int L123 = L123_min; L123 <= L123_max; L123++) {
            double part3_t = 0.0;
            int Tri_x = Tri(2 * L123, 2 * jaco.L_0, 2 * L) +
                        Tri(2 * L123, S_2, jaco.J_2);
            if (Tri_x <= 0)
              continue;
            // cout<<"check point 2"<<endl;
            part3_t = sixJ.quick(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2, S_2,
                                 jaco.J_2);
            // cout<<" C = "<<part3_t<<endl;
            if (part3_t == 0.0)
              continue;
            // cout<<"check point 3"<<endl;
            part3_t =
                part3_t * ninJ.quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123, Sab_2,
                                     1, S_2, jaco.J12_2, jaco.j3_2, jaco.J_2);
            // cout<<" C*D = "<<part3_t<<endl;
            if (part3_t == 0.0)
              continue;
            part3_t = part3_t * (2 * L123 + 1);
            if (part3_t != 0.0) {
              // cout<<"check point 2"<<endl;
            }
            // part4
            double part4 = 0.0;
            for (int L1 = L1_min; L1 <= L1_max; L1++) {
              double part4_t = 0.0;
              if (phase(L1 + jaco.l12) != phase(l_a + l_b))
                continue;
              int N1_2 = Eab - L1 - E12;
              if (N1_2 < 0 || phase(N1_2) < 0)
                continue;
              int N1 = N1_2 / 2;
              int E1 = 2 * N1 + L1;
              part4_t = hobra_x.cal(E1, L1, E12, jaco.l12, E_a, l_a, E_b, l_b,
                                    Lab, d1);
              // cout<<"part4 = "<<part4_t<<endl;
              if (part4_t == 0)
                continue;
              // part5
              double part5 = 0.0;
              for (int L1c = L1c_min; L1c <= L1c_max; L1c++) {
                double part5_t = 0.0;
                if (phase(jaco.L_0 + jaco.l3) != phase(L1 + l_c))
                  continue;
                if ((E0 + E3) != (E1 + E_c))
                  continue;
                int Tri_x = Tri(2 * L1, 2 * l_c, 2 * L1c) +
                            Tri(2 * L, 2 * jaco.l12, 2 * L1c);
                if (Tri_x <= 0)
                  continue;
                part5_t = (2 * L1c + 1) * phase(L1c) *
                          sixJ.quick(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L,
                                     2 * l_c, 2 * L1c);
                if (part5_t == 0.0)
                  continue;
                part5_t =
                    part5_t * sixJ.quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
                                         2 * jaco.L_0, 2 * L, 2 * L1c);
                if (part5_t == 0.0)
                  continue;
                // cout<<"part5_t no hobra :"<<part5_t<<endl;
                part5_t = part5_t * hobra_x.cal(E0, jaco.L_0, E3, jaco.l3, E1,
                                                L1, E_c, l_c, L1c, d2);
                // cout<<"part5_t :"<<part5_t<<endl;
                part5 += part5_t;
              } // part5
              // cout<<"part5"<<part5<<endl;
              part4 += part4_t * part5;
            } // part4
            part3 += part3_t * part4;
          } // part3
          part2 += part2_t * part3;
        } // L
      }   // Lab
    }     // S_2
  }       // Sab_2

  // cout<<part1<<"\t"<<part2<<endl;
  result = part1 * part2;
  return result;
}

double TransCoeff::cal2(State_Jacobi &jaco, State_TP &triP) {
  double result = 0.0;
  if (jaco.T_2 != triP.T_2 || jaco.J2_NL != triP.J_2)
    return 0.0;
  if (jaco.E_NL != triP.E)
    return 0.0;
  if (jaco.T12_2 != triP.Tab_2)
    return 0.0;
  // cout<<"good"<<endl;

  int a_x = triP.a;
  int b_x = triP.b;
  int c_x = triP.c;
  int l_a = SP.state[a_x].l;
  int l_b = SP.state[b_x].l;
  int l_c = SP.state[c_x].l;
  int E_a = SP.state[a_x].E;
  int E_b = SP.state[b_x].E;
  int E_c = SP.state[c_x].E;
  int Eab = SP.state[a_x].E + SP.state[b_x].E;
  int E12 = 2 * jaco.n12 + jaco.l12;
  int E3 = 2 * jaco.n3 + jaco.l3;
  int E0 = 2 * jaco.N_0 + jaco.L_0;

  // cout<<"a :"<<a_x<<endl;
  // cout<<"b :"<<b_x<<endl;
  // cout<<"c :"<<c_x<<endl;

  double d1 = 1.0;
  double d2 = 2.0;

  int L1c_min = std::abs(jaco.L_0 - jaco.l3);
  int L1c_max = jaco.L_0 + jaco.l3;

  int L123_min = std::abs(jaco.l12 - jaco.l3);
  int L123_max = jaco.l12 + jaco.l3;

  // part1
  double part1 = (SP.state[a_x].j_2 + 1) * (SP.state[b_x].j_2 + 1) *
                 (SP.state[c_x].j_2 + 1) * (triP.Jab_2 + 1);
  part1 = part1 * (jaco.J_2 + 1) * (jaco.J12_2 + 1) * (jaco.j3_2 + 1);
  part1 = std::sqrt(part1);
  int phase_x = SP.state[c_x].l + jaco.l12 + jaco.L_0;
  part1 = part1 * phase(phase_x);
  // cout<<"part1 = "<<part1<<endl;
  // part2 part3 ...
  int Lab_min = std::abs(l_a - l_b);
  int Lab_max = l_a + l_b;

  double part2 = 0.0;
  for (int Lab = Lab_min; Lab <= Lab_max; Lab++) {
    int L_min = std::abs(Lab - l_c);
    int L_max = Lab + l_c;
    // cout<<"L : "<<L_min<<"\t"<<L_max<<endl;
    for (int L = L_min; L <= L_max; L++) {
      double part2_t = 0.0;
      int phase_t = Lab + L;
      part2_t = (2 * Lab + 1) * (2 * L + 1) * phase(phase_t);
      // part3 ~~
      int L1_min = std::abs(Lab - jaco.l12);
      int L1_max = Lab + jaco.l12;
      double part3 = 0.0;
      for (int L1 = L1_min; L1 <= L1_max; L1++) {
        double part3_t = 0.0;
        int N1_2 = Eab - L1 - E12;
        if (N1_2 < 0 || phase(N1_2) < 0)
          continue;
        int N1 = N1_2 / 2;
        int E1 = 2 * N1 + L1;
        part3_t =
            hobra_x.cal(E1, L1, E12, jaco.l12, E_a, l_a, E_b, l_b, Lab, d1);
        // part4~~
        double part4 = 0.0;
        for (int L1c = L1c_min; L1c <= L1c_max; L1c++) {
          double part4_t = 0.0;
          part4_t = sixJ.quick(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L, 2 * l_c,
                               2 * L1c);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * hobra_x.cal(E0, jaco.L_0, E3, jaco.l3, E1, L1,
                                          E_c, l_c, L1c, d2);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * (2 * L1c + 1) * phase(L1c);
          // part5
          double part5 = 0.0;
          for (int L123 = L123_min; L123 <= L123_max; L123++) {
            double part5_t = 0.0;
            part5_t = sixJ.quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
                                 2 * jaco.L_0, 2 * L, 2 * L1c);
            if (part5_t == 0)
              continue;
            part5_t = part5_t * (2 * L123 + 1);
            // part6
            double part6 = 0.0;
            for (int Sab_2 = 0; Sab_2 <= 2; Sab_2 = Sab_2 + 2) {
              double part6_t = 0.0;
              part6_t =
                  ninJ.quick(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, Sab_2,
                             SP.state[a_x].j_2, SP.state[b_x].j_2, triP.Jab_2);
              if (part5_t == 0)
                continue;
              part6_t = part6_t * std::sqrt(Sab_2 + 1);
              // part7

              double part7 = 0.0;
              int S2_min = std::abs(Sab_2 - 1);
              int S2_max = Sab_2 + 1;
              for (int S_2 = S2_min; S_2 <= S2_max; S_2 = S_2 + 2) {
                double part7_t = 0.0;
                part7_t = sixJ.quick(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2,
                                     S_2, jaco.J_2);
                if (part7_t == 0)
                  continue;
                part7_t = part7_t * ninJ.quick(2 * Lab, 2 * l_c, 2 * L, Sab_2,
                                               1, S_2, triP.Jab_2,
                                               SP.state[c_x].j_2, triP.J_2);
                if (part7_t == 0)
                  continue;
                part7_t = part7_t * ninJ.quick(2 * jaco.l12, 2 * jaco.l3,
                                               2 * L123, Sab_2, 1, S_2,
                                               jaco.J12_2, jaco.j3_2, jaco.J_2);
                if (part7_t == 0)
                  continue;
                phase_t = (S_2 + jaco.J_2) / 2;
                part7_t = part7_t * (S_2 + 1) * phase(phase_t);

                part7 += part7_t;
              }
              part6 += part6_t * part7;
            }
            part5 += part5_t * part6;
          }
          part4 += part4_t * part5;
        }
        part3 += part3_t * part4;
      }
      part2 += part2_t * part3;
    }
  }
  result = part1 * part2;
  return result;
}

double TransCoeff::cal3(State_Jacobi &jaco, State_TP &triP) {
  // change the phase along with Toku
  double result = 0.0;
  if (jaco.T_2 != triP.T_2 || jaco.J2_NL != triP.J_2)
    return 0.0;
  if (jaco.E_NL != triP.E)
    return 0.0;
  if (jaco.T12_2 != triP.Tab_2)
    return 0.0;
  // cout<<"good"<<endl;

  int a_x = triP.a;
  int b_x = triP.b;
  int c_x = triP.c;

  // cout<<"a :"<<a_x<<endl;
  // cout<<"b :"<<b_x<<endl;
  // cout<<"c :"<<c_x<<endl;

  int l_a = SP.state[a_x].l;
  int l_b = SP.state[b_x].l;
  int l_c = SP.state[c_x].l;
  int E_a = SP.state[a_x].E;
  int E_b = SP.state[b_x].E;
  int E_c = SP.state[c_x].E;
  int Eab = SP.state[a_x].E + SP.state[b_x].E;
  int E12 = 2 * jaco.n12 + jaco.l12;
  int E3 = 2 * jaco.n3 + jaco.l3;
  int E0 = 2 * jaco.N_0 + jaco.L_0;

  double d1 = 1.0;
  double d2 = 2.0;

  int L1c_min = std::abs(jaco.L_0 - jaco.l3);
  int L1c_max = jaco.L_0 + jaco.l3;

  int L123_min = std::abs(jaco.l12 - jaco.l3);
  int L123_max = jaco.l12 + jaco.l3;

  // part1
  double part1 = (SP.state[a_x].j_2 + 1) * (SP.state[b_x].j_2 + 1) *
                 (SP.state[c_x].j_2 + 1) * (triP.Jab_2 + 1);
  part1 = part1 * (jaco.J_2 + 1) * (jaco.J12_2 + 1) * (jaco.j3_2 + 1);
  part1 = std::sqrt(part1);
  //  int phase_x = SP.state[c_x].l+jaco.l12+jaco.L_0;
  int phase_x = SP.state[c_x].l + jaco.l12;
  part1 = part1 * phase(phase_x);
  // cout<<"part1 = "<<part1<<endl;
  // part2 part3 ...
  int Lab_min = std::abs(l_a - l_b);
  int Lab_max = l_a + l_b;

  double part2 = 0.0;
  for (int Lab = Lab_min; Lab <= Lab_max; Lab++) {
    int L_min = std::abs(Lab - l_c);
    int L_max = Lab + l_c;
    // cout<<"L : "<<L_min<<"\t"<<L_max<<endl;
    for (int L = L_min; L <= L_max; L++) {
      double part2_t = 0.0;
      int phase_t = Lab + L;
      part2_t = (2 * Lab + 1) * (2 * L + 1) * phase(phase_t);
      // part3 ~~
      int L1_min = std::abs(Lab - jaco.l12);
      int L1_max = Lab + jaco.l12;
      double part3 = 0.0;
      for (int L1 = L1_min; L1 <= L1_max; L1++) {
        double part3_t = 0.0;
        int N1_2 = Eab - L1 - E12;
        if (N1_2 < 0 || phase(N1_2) < 0)
          continue;
        int N1 = N1_2 / 2;
        int E1 = 2 * N1 + L1;
        part3_t =
            hobra_x.cal(E1, L1, E12, jaco.l12, E_a, l_a, E_b, l_b, Lab, d1);
        // part4~~
        double part4 = 0.0;
        for (int L1c = L1c_min; L1c <= L1c_max; L1c++) {
          double part4_t = 0.0;
          part4_t = sixJ.quick(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L, 2 * l_c,
                               2 * L1c);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * hobra_x.cal(E0, jaco.L_0, E3, jaco.l3, E1, L1,
                                          E_c, l_c, L1c, d2);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * (2 * L1c + 1) * phase(L1c);
          // part5
          double part5 = 0.0;
          for (int L123 = L123_min; L123 <= L123_max; L123++) {
            double part5_t = 0.0;
            part5_t = sixJ.quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
                                 2 * jaco.L_0, 2 * L, 2 * L1c);
            if (part5_t == 0)
              continue;
            part5_t = part5_t * (2 * L123 + 1);
            // part6
            double part6 = 0.0;
            for (int Sab_2 = 0; Sab_2 <= 2; Sab_2 = Sab_2 + 2) {
              double part6_t = 0.0;
              part6_t =
                  ninJ.quick(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, Sab_2,
                             SP.state[a_x].j_2, SP.state[b_x].j_2, triP.Jab_2);
              if (part5_t == 0)
                continue;
              part6_t = part6_t * std::sqrt(Sab_2 + 1);
              // part7

              double part7 = 0.0;
              int S2_min = std::abs(Sab_2 - 1);
              int S2_max = Sab_2 + 1;
              for (int S_2 = S2_min; S_2 <= S2_max; S_2 = S_2 + 2) {
                double part7_t = 0.0;
                part7_t = sixJ.quick(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2,
                                     S_2, jaco.J_2);
                if (part7_t == 0)
                  continue;
                part7_t = part7_t * ninJ.quick(2 * Lab, 2 * l_c, 2 * L, Sab_2,
                                               1, S_2, triP.Jab_2,
                                               SP.state[c_x].j_2, triP.J_2);
                if (part7_t == 0)
                  continue;
                part7_t = part7_t * ninJ.quick(2 * jaco.l12, 2 * jaco.l3,
                                               2 * L123, Sab_2, 1, S_2,
                                               jaco.J12_2, jaco.j3_2, jaco.J_2);
                if (part7_t == 0)
                  continue;
                phase_t = (S_2 + triP.J_2) / 2;
                part7_t = part7_t * (S_2 + 1) * phase(phase_t);

                part7 += part7_t;
              }
              part6 += part6_t * part7;
            }
            part5 += part5_t * part6;
          }
          part4 += part4_t * part5;
        }
        part3 += part3_t * part4;
      }
      part2 += part2_t * part3;
    }
  }
  result = part1 * part2;
  return result;
}

double TransCoeff::cal4(State_Jacobi &jaco, State_TP &triP) {
  // change the phase along with Toku
  double result = 0.0;
  if (jaco.T_2 != triP.T_2 || jaco.J2_NL != triP.J_2)
    return 0.0;
  if (jaco.E_NL != triP.E)
    return 0.0;
  if (jaco.T12_2 != triP.Tab_2)
    return 0.0;
  // cout<<"good"<<endl;

  int a_x = triP.a;
  int b_x = triP.b;
  int c_x = triP.c;
  int l_a = SP.state[a_x].l;
  int l_b = SP.state[b_x].l;
  int l_c = SP.state[c_x].l;
  int E_a = SP.state[a_x].E;
  int E_b = SP.state[b_x].E;
  int E_c = SP.state[c_x].E;
  int Eab = SP.state[a_x].E + SP.state[b_x].E;
  int E12 = 2 * jaco.n12 + jaco.l12;
  int E3 = 2 * jaco.n3 + jaco.l3;
  int E0 = 2 * jaco.N_0 + jaco.L_0;

  // cout<<"a :"<<a_x<<endl;
  // cout<<"b :"<<b_x<<endl;
  // cout<<"c :"<<c_x<<endl;

  double d1 = 1.0;
  double d2 = 2.0;

  int L1c_min = std::abs(jaco.L_0 - jaco.l3);
  int L1c_max = jaco.L_0 + jaco.l3;

  int L123_min = std::abs(jaco.l12 - jaco.l3);
  int L123_max = jaco.l12 + jaco.l3;

  // part1
  double part1 = (SP.state[a_x].j_2 + 1) * (SP.state[b_x].j_2 + 1) *
                 (SP.state[c_x].j_2 + 1) * (triP.Jab_2 + 1);
  part1 = part1 * (jaco.J_2 + 1) * (jaco.J12_2 + 1) * (jaco.j3_2 + 1);
  part1 = std::sqrt(part1);
  //  int phase_x = SP.state[c_x].l+jaco.l12+jaco.L_0;
  int phase_x = SP.state[c_x].l + jaco.l12;
  part1 = part1 * phase(phase_x);
  // cout<<"part1 = "<<part1<<endl;
  // part2 part3 ...
  int Lab_min = std::abs(l_a - l_b);
  int Lab_max = l_a + l_b;

  double part2 = 0.0;
  for (int Lab = Lab_min; Lab <= Lab_max; Lab++) {
    int L_min = std::abs(Lab - l_c);
    int L_max = Lab + l_c;
    // cout<<"L : "<<L_min<<"\t"<<L_max<<endl;
    for (int L = L_min; L <= L_max; L++) {
      double part2_t = 0.0;
      int phase_t = Lab + L;
      part2_t = (2 * Lab + 1) * (2 * L + 1) * phase(phase_t);
      // part3 ~~
      int L1_min = std::abs(Lab - jaco.l12);
      int L1_max = Lab + jaco.l12;
      double part3 = 0.0;
      for (int L1 = L1_min; L1 <= L1_max; L1++) {
        double part3_t = 0.0;
        int N1_2 = Eab - L1 - E12;
        if (N1_2 < 0 || phase(N1_2) < 0)
          continue;
        int N1 = N1_2 / 2;
        int E1 = 2 * N1 + L1;
        part3_t =
            hobra_x.cal(E1, L1, E12, jaco.l12, E_a, l_a, E_b, l_b, Lab, d1);
        // part4~~
        double part4 = 0.0;
        for (int L1c = L1c_min; L1c <= L1c_max; L1c++) {
          double part4_t = 0.0;
          part4_t = sixJ.quick(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L, 2 * l_c,
                               2 * L1c);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * hobra_x.cal(E0, jaco.L_0, E3, jaco.l3, E1, L1,
                                          E_c, l_c, L1c, d2);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * (2 * L1c + 1) * phase(L1c);
          // part5
          double part5 = 0.0;
          for (int L123 = L123_min; L123 <= L123_max; L123++) {
            double part5_t = 0.0;
            part5_t = sixJ.quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
                                 2 * jaco.L_0, 2 * L, 2 * L1c);
            if (part5_t == 0)
              continue;
            part5_t = part5_t * (2 * L123 + 1);
            // part6
            double part6 = 0.0;
            for (int Sab_2 = 0; Sab_2 <= 2; Sab_2 = Sab_2 + 2) {
              double part6_t = 0.0;
              part6_t =
                  ninJ.quick(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, Sab_2,
                             SP.state[a_x].j_2, SP.state[b_x].j_2, triP.Jab_2);
              if (part5_t == 0)
                continue;
              part6_t = part6_t * std::sqrt(Sab_2 + 1);
              // part7

              double part7 = 0.0;
              int S2_min = std::abs(Sab_2 - 1);
              int S2_max = Sab_2 + 1;
              for (int S_2 = S2_min; S_2 <= S2_max; S_2 = S_2 + 2) {
                double part7_t = 0.0;
                part7_t = sixJ.quick(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2,
                                     S_2, jaco.J_2);
                if (part7_t == 0)
                  continue;
                part7_t = part7_t * ninJ.quick(2 * Lab, 2 * l_c, 2 * L, Sab_2,
                                               1, S_2, triP.Jab_2,
                                               SP.state[c_x].j_2, triP.J_2);
                if (part7_t == 0)
                  continue;
                part7_t = part7_t * ninJ.quick(2 * jaco.l12, 2 * jaco.l3,
                                               2 * L123, Sab_2, 1, S_2,
                                               jaco.J12_2, jaco.j3_2, jaco.J_2);
                if (part7_t == 0)
                  continue;
                phase_t = (S_2 + triP.J_2) / 2;
                part7_t = part7_t * (S_2 + 1) * phase(phase_t);

                part7 += part7_t;
              }
              part6 += part6_t * part7;
            }
            part5 += part5_t * part6;
          }
          part4 += part4_t * part5;
        }
        part3 += part3_t * part4;
      }
      part2 += part2_t * part3;
    }
  }
  result = part1 * part2;
  return result;
}
double TransCoeff::cal5(State_Jacobi &jaco, State_TP &triP) {
  // change the phase along with Toku
  double result = 0.0;
  if (jaco.T_2 != triP.T_2 || jaco.J2_NL != triP.J_2)
    return 0.0;
  if (jaco.E_NL != triP.E)
    return 0.0;
  if (jaco.T12_2 != triP.Tab_2)
    return 0.0;
  // cout<<"good"<<endl;

  int a_x = triP.a;
  int b_x = triP.b;
  int c_x = triP.c;

  // cout<<"a :"<<a_x<<endl;
  // cout<<"b :"<<b_x<<endl;
  // cout<<"c :"<<c_x<<endl;

  int l_a = SP.state[a_x].l;
  int l_b = SP.state[b_x].l;
  int l_c = SP.state[c_x].l;
  int E_a = SP.state[a_x].E;
  int E_b = SP.state[b_x].E;
  int E_c = SP.state[c_x].E;
  int Eab = SP.state[a_x].E + SP.state[b_x].E;
  int E12 = 2 * jaco.n12 + jaco.l12;
  int E3 = 2 * jaco.n3 + jaco.l3;
  int E0 = 2 * jaco.N_0 + jaco.L_0;
  int S12_2 = jaco.S12_2;

  double d1 = 1.0;
  double d2 = 2.0;

  int L1c_min = std::abs(jaco.L_0 - jaco.l3);
  int L1c_max = jaco.L_0 + jaco.l3;

  int L123_min = std::abs(jaco.l12 - jaco.l3);
  int L123_max = jaco.l12 + jaco.l3;

  int S2_min = std::abs(S12_2 - 1);
  int S2_max = S12_2 + 1;

  // part1
  double part1 = (SP.state[a_x].j_2 + 1) * (SP.state[b_x].j_2 + 1) *
                 (SP.state[c_x].j_2 + 1) * (triP.Jab_2 + 1);
  part1 =
      part1 * (jaco.J_2 + 1) * (S12_2 + 1) * (jaco.J12_2 + 1) * (jaco.j3_2 + 1);
  part1 = std::sqrt(part1);
  //  int phase_x = SP.state[c_x].l+jaco.l12+jaco.L_0;
  int phase_x = SP.state[c_x].l + jaco.l12;
  part1 = part1 * phase(phase_x);
  // cout<<"part1 = "<<part1<<endl;
  // part2 part3 ...
  int Lab_min = std::abs(l_a - l_b);
  int Lab_max = l_a + l_b;

  double part2 = 0.0;
  for (int Lab = Lab_min; Lab <= Lab_max; Lab++) { // range is Small ~ 0-10
    int L_min = std::abs(Lab - l_c);
    int L_max = Lab + l_c;
    // cout<<"L : "<<L_min<<"\t"<<L_max<<endl;
    for (int L = L_min; L <= L_max; L++) { // range is Small ~ 0-10
      double part2_t = 0.0;
      int phase_t = Lab + L;
      part2_t = ninJ.quick(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, S12_2,
                           SP.state[a_x].j_2, SP.state[b_x].j_2, triP.Jab_2);
      if (part2_t == 0.0)
        continue;
      part2_t = part2_t * (2 * Lab + 1) * (2 * L + 1) * phase(phase_t);
      // part3 ~~
      int L1_min = std::abs(Lab - jaco.l12);
      int L1_max = Lab + jaco.l12;
      double part3 = 0.0;
      for (int L1 = L1_min; L1 <= L1_max; L1++) { // range is Small ~ 0-10
        double part3_t = 0.0;
        int N1_2 = Eab - L1 - E12;
        if (N1_2 < 0 || phase(N1_2) < 0)
          continue;
        int N1 = N1_2 / 2;
        int E1 = 2 * N1 + L1;
        part3_t =
            hobra_x.cal(E1, L1, E12, jaco.l12, E_a, l_a, E_b, l_b, Lab, d1);
        // part4~~
        double part4 = 0.0;
        for (int L1c = L1c_min; L1c <= L1c_max;
             L1c++) { // range is Small ~ 0-10
          double part4_t = 0.0;
          part4_t = sixJ.quick(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L, 2 * l_c,
                               2 * L1c);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * hobra_x.cal(E0, jaco.L_0, E3, jaco.l3, E1, L1,
                                          E_c, l_c, L1c, d2);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * (2 * L1c + 1) * phase(L1c);
          // part5
          double part5 = 0.0;
          for (int L123 = L123_min; L123 <= L123_max;
               L123++) { // range is Small ~ 0-10
            double part5_t = 0.0;
            part5_t = sixJ.quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
                                 2 * jaco.L_0, 2 * L, 2 * L1c);
            if (part5_t == 0)
              continue;
            part5_t = part5_t * (2 * L123 + 1);
            // part6
            double part6 = 0.0;
            for (int S_2 = S2_min; S_2 <= S2_max;
                 S_2 = S_2 + 2) { // range is Small ~ 0-5
              double part6_t = 0.0;
              part6_t = sixJ.quick(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2, S_2,
                                   jaco.J_2);
              if (part6_t == 0)
                continue;
              part6_t =
                  part6_t * ninJ.quick(2 * Lab, 2 * l_c, 2 * L, S12_2, 1, S_2,
                                       triP.Jab_2, SP.state[c_x].j_2, triP.J_2);
              if (part6_t == 0)
                continue;
              part6_t = part6_t * ninJ.quick(2 * jaco.l12, 2 * jaco.l3,
                                             2 * L123, S12_2, 1, S_2,
                                             jaco.J12_2, jaco.j3_2, jaco.J_2);
              if (part6_t == 0)
                continue;
              phase_t = (S_2 + triP.J_2) / 2;
              part6_t = part6_t * (S_2 + 1) * phase(phase_t);

              part6 += part6_t;
            }
            part5 += part5_t * part6;
          }
          part4 += part4_t * part5;
        }
        part3 += part3_t * part4;
      }
      part2 += part2_t * part3;
    }
  }
  result = part1 * part2;
  return result;
}

//#pragma acc routine(sqrt) seq

// using std::sqrt;
// using std::fabs;
//
// #pragma acc routine seq
// template <typename T> T sqabs(T x) { return sqrt(fabs(x)); }
//#pragma acc routine vector
double TransCoeff::cal5_OA(State_Jacobi &jaco, State_TP &triP) {
  // change the phase along with Toku
  double result = 0.0;
  // return result;
  if (jaco.T_2 != triP.T_2 || jaco.J2_NL != triP.J_2)
    return 0.0;
  if (jaco.E_NL != triP.E)
    return 0.0;
  if (jaco.T12_2 != triP.Tab_2)
    return 0.0;
  // cout<<"good"<<endl;

  int a_x = triP.a;
  int b_x = triP.b;
  int c_x = triP.c;

  // cout<<"a :"<<a_x<<endl;
  // cout<<"b :"<<b_x<<endl;
  // cout<<"c :"<<c_x<<endl;

  int l_a = SP.state[a_x].l;
  int l_b = SP.state[b_x].l;
  int l_c = SP.state[c_x].l;
  int E_a = SP.state[a_x].E;
  int E_b = SP.state[b_x].E;
  int E_c = SP.state[c_x].E;
  int Eab = SP.state[a_x].E + SP.state[b_x].E;
  int E12 = 2 * jaco.n12 + jaco.l12;
  int E3 = 2 * jaco.n3 + jaco.l3;
  int E0 = 2 * jaco.N_0 + jaco.L_0;
  int S12_2 = jaco.S12_2;

  double d1 = 1.0;
  double d2 = 2.0;

  int L1c_min = std::abs(jaco.L_0 - jaco.l3);
  int L1c_max = jaco.L_0 + jaco.l3;

  int L123_min = std::abs(jaco.l12 - jaco.l3);
  int L123_max = jaco.l12 + jaco.l3;

  int S2_min = std::abs(S12_2 - 1);
  int S2_max = S12_2 + 1;

  // part1
  double part1 = (SP.state[a_x].j_2 + 1) * (SP.state[b_x].j_2 + 1) *
                 (SP.state[c_x].j_2 + 1) * (triP.Jab_2 + 1);
  part1 =
      part1 * (jaco.J_2 + 1) * (S12_2 + 1) * (jaco.J12_2 + 1) * (jaco.j3_2 + 1);

  part1 = std::sqrt(part1);
  // part1 = cal_sqrt(part1);

  //  int phase_x = SP.state[c_x].l+jaco.l12+jaco.L_0;
  int phase_x = SP.state[c_x].l + jaco.l12;
  part1 = part1 * phase(phase_x);
  // cout<<"part1 = "<<part1<<endl;
  // part2 part3 ...
  int Lab_min = std::abs(l_a - l_b);
  int Lab_max = l_a + l_b;

  double part2 = 0.0;
  //#pragma acc loop vector reduction(+ : part2)
  nvtxRangePush("Lab loop");
  for (int Lab = Lab_min; Lab <= Lab_max; Lab++) { // range is Small ~ 0-10
    int L_min = std::abs(Lab - l_c);
    int L_max = Lab + l_c;
    // cout<<"L : "<<L_min<<"\t"<<L_max<<endl;
    nvtxRangePush("L loop");
    for (int L = L_min; L <= L_max; L++) { // range is Small ~ 0-10
      double part2_t = 0.0;
      int phase_t = Lab + L;
      part2_t = ninJ.quick(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, S12_2,
                           SP.state[a_x].j_2, SP.state[b_x].j_2, triP.Jab_2);
      // part2_t = NineJSymbol_2(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, S12_2,
      //                         SP.state[a_x].j_2, SP.state[b_x].j_2,
      //                         triP.Jab_2);
      // part2_t = ninJ.cal(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, S12_2,
      //                     SP.state[a_x].j_2, SP.state[b_x].j_2, triP.Jab_2);
      if (part2_t == 0.0)
        continue;
      part2_t = part2_t * (2 * Lab + 1) * (2 * L + 1) * phase(phase_t);
      // part3 ~~
      int L1_min = std::abs(Lab - jaco.l12);
      int L1_max = Lab + jaco.l12;
      double part3 = 0.0;
      nvtxRangePush("L1 loop");
      for (int L1 = L1_min; L1 <= L1_max; L1++) { // range is Small ~ 0-10
        double part3_t = 0.0;
        int N1_2 = Eab - L1 - E12;
        if (N1_2 < 0 || phase(N1_2) < 0)
          continue;
        int N1 = N1_2 / 2;
        int E1 = 2 * N1 + L1;
        part3_t =
            hobra_x.cal(E1, L1, E12, jaco.l12, E_a, l_a, E_b, l_b, Lab, d1);
        // part4~~
        double part4 = 0.0;
        nvtxRangePush("L1c loop");
        for (int L1c = L1c_min; L1c <= L1c_max;
             L1c++) { // range is Small ~ 0-10
          double part4_t = 0.0;
          part4_t = sixJ.quick(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L, 2 * l_c,
                               2 * L1c);
          // part4_t = SixJSymbol_2(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L, 2 *
          // l_c,
          //                        2 * L1c);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * hobra_x.cal(E0, jaco.L_0, E3, jaco.l3, E1, L1,
                                          E_c, l_c, L1c, d2);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * (2 * L1c + 1) * phase(L1c);
          // part5
          double part5 = 0.0;
          nvtxRangePush("L123 loop");
          for (int L123 = L123_min; L123 <= L123_max;
               L123++) { // range is Small ~ 0-10
            double part5_t = 0.0;
            part5_t = sixJ.quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
                                 2 * jaco.L_0, 2 * L, 2 * L1c);
            // part5_t = SixJSymbol_2(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
            //                        2 * jaco.L_0, 2 * L, 2 * L1c);
            if (part5_t == 0)
              continue;
            part5_t = part5_t * (2 * L123 + 1);
            // part6
            double part6 = 0.0;
            for (int S_2 = S2_min; S_2 <= S2_max;
                 S_2 = S_2 + 2) { // range is Small ~ 0-5
              double part6_t = 0.0;
              part6_t = sixJ.quick(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2, S_2,
                                   jaco.J_2);
              // part6_t = SixJSymbol_2(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2,
              //                        S_2, jaco.J_2);
              if (part6_t == 0)
                continue;
              part6_t =
                  part6_t * ninJ.quick(2 * Lab, 2 * l_c, 2 * L, S12_2, 1, S_2,
                                       triP.Jab_2, SP.state[c_x].j_2, triP.J_2);
              // part6_t = part6_t * NineJSymbol_2(2 * Lab, 2 * l_c, 2 * L,
              // S12_2,
              //                                   1, S_2, triP.Jab_2,
              //                                   SP.state[c_x].j_2, triP.J_2);
              if (part6_t == 0)
                continue;
              part6_t = part6_t * ninJ.quick(2 * jaco.l12, 2 * jaco.l3,
                                             2 * L123, S12_2, 1, S_2,
                                             jaco.J12_2, jaco.j3_2, jaco.J_2);
              // part6_t = part6_t * NineJSymbol_2(2 * jaco.l12, 2 * jaco.l3,
              //                                2 * L123, S12_2, 1, S_2,
              //                                jaco.J12_2, jaco.j3_2,
              //                                jaco.J_2);
              if (part6_t == 0)
                continue;
              phase_t = (S_2 + triP.J_2) / 2;
              part6_t = part6_t * (S_2 + 1) * phase(phase_t);

              part6 += part6_t;
            }
            part5 += part5_t * part6;
          } // L123
          nvtxRangePop();
          part4 += part4_t * part5;
        } // L1c
        nvtxRangePop();
        part3 += part3_t * part4;
      } // L1
      nvtxRangePop();
      part2 += part2_t * part3;
    } //   L
    nvtxRangePop();
  } // Lab
  nvtxRangePop();
  result = part1 * part2;
  return result;
}
//#pragma acc routine vector
double TransCoeff::cal5_OA(State_Jacobi &jaco, State_TP &triP,
                           vector<State_SP> &SP_state) {
  // change the phase along with Toku
  double result = 0.0;
  //
  // return result;
  //
  if (jaco.T_2 != triP.T_2 || jaco.J2_NL != triP.J_2)
    return 0.0;
  if (jaco.E_NL != triP.E)
    return 0.0;
  if (jaco.T12_2 != triP.Tab_2)
    return 0.0;
  // cout<<"good"<<endl;

  int a_x = triP.a;
  int b_x = triP.b;
  int c_x = triP.c;

  // cout<<"a :"<<a_x<<endl;
  // cout<<"b :"<<b_x<<endl;
  // cout<<"c :"<<c_x<<endl;
  int l_a = SP_state[a_x].l;
  int l_b = SP_state[b_x].l;
  int l_c = SP_state[c_x].l;
  int E_a = SP_state[a_x].E;
  int E_b = SP_state[b_x].E;
  int E_c = SP_state[c_x].E;
  int Eab = SP_state[a_x].E + SP_state[b_x].E;
  int E12 = 2 * jaco.n12 + jaco.l12;
  int E3 = 2 * jaco.n3 + jaco.l3;
  int E0 = 2 * jaco.N_0 + jaco.L_0;
  int S12_2 = jaco.S12_2;

  double d1 = 1.0;
  double d2 = 2.0;

  int L1c_min = std::abs(jaco.L_0 - jaco.l3);
  int L1c_max = jaco.L_0 + jaco.l3;

  int L123_min = std::abs(jaco.l12 - jaco.l3);
  int L123_max = jaco.l12 + jaco.l3;

  int S2_min = std::abs(S12_2 - 1);
  int S2_max = S12_2 + 1;

  // return result;

  // part1
  double part1 = (SP_state[a_x].j_2 + 1) * (SP_state[b_x].j_2 + 1) *
                 (SP_state[c_x].j_2 + 1) * (triP.Jab_2 + 1);
  part1 =
      part1 * (jaco.J_2 + 1) * (S12_2 + 1) * (jaco.J12_2 + 1) * (jaco.j3_2 + 1);

  part1 = std::sqrt(part1);
  // part1 = cal_sqrt(part1);

  //  int phase_x = SP.state[c_x].l+jaco.l12+jaco.L_0;
  int phase_x = SP_state[c_x].l + jaco.l12;
  part1 = part1 * phase(phase_x);
  // cout<<"part1 = "<<part1<<endl;
  // part2 part3 ...
  int Lab_min = std::abs(l_a - l_b);
  int Lab_max = l_a + l_b;

  // return result;
  double part2 = 0.0;
  //#pragma acc loop reduction(+ : part2)
  //#pragma acc parallel loop worker collapse(3)
  for (int Lab = Lab_min; Lab <= Lab_max; Lab++) { // range is Small ~ 0-10
    int L_min = std::abs(Lab - l_c);
    int L_max = Lab + l_c;
    //#pragma acc loop reduction(+ : part2)
    for (int L = L_min; L <= L_max; L++) { // range is Small ~ 0-10
      double part2_t = 0.0;
      int phase_t = Lab + L;

      part2_t = 1; // ninJ.quick(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, S12_2,
                   //       SP_state[a_x].j_2, SP_state[b_x].j_2, triP.Jab_2);

      if (part2_t == 0.0)
        continue;
      part2_t = part2_t * (2 * Lab + 1) * (2 * L + 1) * phase(phase_t);
      // part3 ~~
      int L1_min = std::abs(Lab - jaco.l12);
      int L1_max = Lab + jaco.l12;
      double part3 = 0.0;
      //#pragma acc loop reduction(+ : part2)
      for (int L1 = L1_min; L1 <= L1_max; L1++) { // range is Small ~ 0-10

        double part3_t = 0.0;
        int N1_2 = Eab - L1 - E12;
        if (N1_2 < 0 || phase(N1_2) < 0)
          continue;
        int N1 = N1_2 / 2;
        int E1 = 2 * N1 + L1;
        part3_t =
            hobra_x.cal2(E1, L1, E12, jaco.l12, E_a, l_a, E_b, l_b, Lab, d1);
        // part4~~
        double part4 = 0.0;
        // #pragma acc loop reduction(+ : part4)
        for (int L1c = L1c_min; L1c <= L1c_max;
             L1c++) { // range is Small ~ 0-10

          double part4_t = 0.0;
          part4_t =
              1; // sixJ.quick(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L, 2 * l_c,
                 //      2 * L1c);
          // part4_t = SixJSymbol_2(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L, 2 *
          // l_c,
          //                        2 * L1c);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * hobra_x.cal2(E0, jaco.L_0, E3, jaco.l3, E1, L1,
                                           E_c, l_c, L1c, d2);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * (2 * L1c + 1) * phase(L1c);
          // part5
          double part5 = 0.0;
          //    nvtxRangePush("L123 loop");
          for (int L123 = L123_min; L123 <= L123_max;
               L123++) { // range is Small ~ 0-10

            double part5_t = 0.0;
            part5_t = 1; // sixJ.quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
                         //       2 * jaco.L_0, 2 * L, 2 * L1c);
            // part5_t = SixJSymbol_2(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
            //                        2 * jaco.L_0, 2 * L, 2 * L1c);
            if (part5_t == 0)
              continue;
            part5_t = part5_t * (2 * L123 + 1);
            // part6
            double part6 = 0.0;
            for (int S_2 = S2_min; S_2 <= S2_max;
                 S_2 = S_2 + 2) { // range is Small ~ 0-5
              // continue;
              double part6_t = 0.0;
              part6_t = 1;
              ; // sixJ.quick(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2, S_2,
                //     jaco.J_2);
              // part6_t = SixJSymbol_2(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2,
              //                        S_2, jaco.J_2);
              if (part6_t == 0)
                continue;
              part6_t = part6_t *
                        1; // ninJ.quick(2 * Lab, 2 * l_c, 2 * L, S12_2, 1, S_2,
                           //       triP.Jab_2, SP_state[c_x].j_2, triP.J_2);
              // part6_t = part6_t * NineJSymbol_2(2 * Lab, 2 * l_c, 2 * L,
              // S12_2,
              //                                   1, S_2, triP.Jab_2,
              //                                   SP.state[c_x].j_2, triP.J_2);
              if (part6_t == 0)
                continue;
              part6_t = part6_t * 1; // ninJ.quick(2 * jaco.l12, 2 * jaco.l3,
                                     //       2 * L123, S12_2, 1, S_2,
                                     //     jaco.J12_2, jaco.j3_2, jaco.J_2);
              // part6_t = part6_t * NineJSymbol_2(2 * jaco.l12, 2 * jaco.l3,
              //                                2 * L123, S12_2, 1, S_2,
              //                                jaco.J12_2, jaco.j3_2,
              //                                jaco.J_2);
              if (part6_t == 0)
                continue;
              phase_t = (S_2 + triP.J_2) / 2;
              part6_t = part6_t * (S_2 + 1) * phase(phase_t);

              part6 += part6_t;
            }
            part5 += part5_t * part6;
          } // nvtxRangePop();
          part4 += part4_t * part5;
        } // nvtxRangePop();
        part3 += part3_t * part4;
      } // nvtxRangePop();
      part2 += part2_t * part3;
    } // nvtxRangePop();
  }   // nvtxRangePop();
  result = part1 * part2;
  return result;
}

#pragma acc routine vector
double TransCoeff::cal5_OA(State_Jacobi &jaco, State_TP &triP,
                           vector<State_SP> &SP_state, Pre_coeff &PC) {
  // change the phase along with Toku
  double result = 0.0;
  //
  // return result;
  //
  if (jaco.T_2 != triP.T_2 || jaco.J2_NL != triP.J_2)
    return 0.0;
  if (jaco.E_NL != triP.E)
    return 0.0;
  if (jaco.T12_2 != triP.Tab_2)
    return 0.0;
  // cout<<"good"<<endl;

  int a_x = triP.a;
  int b_x = triP.b;
  int c_x = triP.c;

  // cout<<"a :"<<a_x<<endl;
  // cout<<"b :"<<b_x<<endl;
  // cout<<"c :"<<c_x<<endl;
  int l_a = SP_state[a_x].l;
  int l_b = SP_state[b_x].l;
  int l_c = SP_state[c_x].l;
  int E_a = SP_state[a_x].E;
  int E_b = SP_state[b_x].E;
  int E_c = SP_state[c_x].E;
  int Eab = SP_state[a_x].E + SP_state[b_x].E;
  int E12 = 2 * jaco.n12 + jaco.l12;
  int E3 = 2 * jaco.n3 + jaco.l3;
  int E0 = 2 * jaco.N_0 + jaco.L_0;
  int S12_2 = jaco.S12_2;

  double d1 = 1.0;
  double d2 = 2.0;

  int L1c_min = std::abs(jaco.L_0 - jaco.l3);
  int L1c_max = jaco.L_0 + jaco.l3;

  int L123_min = std::abs(jaco.l12 - jaco.l3);
  int L123_max = jaco.l12 + jaco.l3;

  int S2_min = std::abs(S12_2 - 1);
  int S2_max = S12_2 + 1;

  // return result;

  // part1
  double part1 = (SP_state[a_x].j_2 + 1) * (SP_state[b_x].j_2 + 1) *
                 (SP_state[c_x].j_2 + 1) * (triP.Jab_2 + 1);
  part1 =
      part1 * (jaco.J_2 + 1) * (S12_2 + 1) * (jaco.J12_2 + 1) * (jaco.j3_2 + 1);

  part1 = std::sqrt(part1);
  int phase_x = SP_state[c_x].l + jaco.l12;
  part1 = part1 * phase(phase_x);
  // cout<<"part1 = "<<part1<<endl;
  // part2 part3 ...
  int Lab_min = std::abs(l_a - l_b);
  int Lab_max = l_a + l_b;

  // return result;
  double part2 = 0.0;
  for (int Lab = Lab_min; Lab <= Lab_max; Lab++) { // range is Small ~ 0-10
    int L_min = std::abs(Lab - l_c);
    int L_max = Lab + l_c;
    //#pragma acc loop reduction(+ : part2)
    for (int L = L_min; L <= L_max; L++) { // range is Small ~ 0-10
      double part2_t = 0.0;
      int phase_t = Lab + L;
      part2_t = PC.NinJ_quick(2 * l_a, 2 * l_b, 2 * Lab, 1, 1, S12_2,
                              SP_state[a_x].j_2, SP_state[b_x].j_2, triP.Jab_2);

      if (part2_t == 0.0)
        continue;
      part2_t = part2_t * (2 * Lab + 1) * (2 * L + 1) * phase(phase_t);
      // part3 ~~
      int L1_min = std::abs(Lab - jaco.l12);
      int L1_max = Lab + jaco.l12;
      double part3 = 0.0;
      for (int L1 = L1_min; L1 <= L1_max; L1++) { // range is Small ~ 0-10

        double part3_t = 0.0;
        int N1_2 = Eab - L1 - E12;
        if (N1_2 < 0 || phase(N1_2) < 0)
          continue;
        int N1 = N1_2 / 2;
        int E1 = 2 * N1 + L1;
        part3_t =
            PC.hobra_cal(E1, L1, E12, jaco.l12, E_a, l_a, E_b, l_b, Lab, d1);
        // part4~~
        double part4 = 0.0;
        // #pragma acc loop reduction(+ : part4)
        for (int L1c = L1c_min; L1c <= L1c_max;
             L1c++) { // range is Small ~ 0-10

          double part4_t = 0.0;
          part4_t = PC.SixJ_quick(2 * L1, 2 * jaco.l12, 2 * Lab, 2 * L, 2 * l_c,
                                  2 * L1c);
          if (part4_t == 0)
            continue;
          // part4_t = part4_t * hobra_x.cal(E0, jaco.L_0, E3, jaco.l3, E1, L1,
          //                                 E_c, l_c, L1c, d2);
          part4_t = part4_t * PC.hobra_cal(E0, jaco.L_0, E3, jaco.l3, E1, L1,
                                           E_c, l_c, L1c, d2);
          if (part4_t == 0)
            continue;
          part4_t = part4_t * (2 * L1c + 1) * phase(L1c);
          // part5
          double part5 = 0.0;
          //    nvtxRangePush("L123 loop");
          for (int L123 = L123_min; L123 <= L123_max;
               L123++) { // range is Small ~ 0-10

            double part5_t = 0.0;
            part5_t = PC.SixJ_quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
                                    2 * jaco.L_0, 2 * L, 2 * L1c);
            // part5_t = sixJ.quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
            //                      2 * jaco.L_0, 2 * L, 2 * L1c);
            // part5_t = SixJSymbol_2(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
            //                        2 * jaco.L_0, 2 * L, 2 * L1c);
            if (part5_t == 0)
              continue;
            part5_t = part5_t * (2 * L123 + 1);
            // part6
            double part6 = 0.0;
            for (int S_2 = S2_min; S_2 <= S2_max;
                 S_2 = S_2 + 2) { // range is Small ~ 0-5
              double part6_t = 0.0;
              part6_t = PC.SixJ_quick(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2,
                                      S_2, jaco.J_2);
              // part6_t = sixJ.quick(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2,
              // S_2,
              //                      jaco.J_2);
              // part6_t = SixJSymbol_2(2 * L123, 2 * jaco.L_0, 2 * L, triP.J_2,
              //                        S_2, jaco.J_2);
              if (part6_t == 0)
                continue;
              part6_t = part6_t * PC.NinJ_quick(2 * Lab, 2 * l_c, 2 * L, S12_2,
                                                1, S_2, triP.Jab_2,
                                                SP_state[c_x].j_2, triP.J_2);
              // part6_t =
              //     part6_t * ninJ.quick(2 * Lab, 2 * l_c, 2 * L, S12_2, 1,
              //     S_2,
              //                          triP.Jab_2, SP_state[c_x].j_2,
              //                          triP.J_2);
              // part6_t = part6_t * NineJSymbol_2(2 * Lab, 2 * l_c, 2 * L,
              // S12_2,
              //                                   1, S_2, triP.Jab_2,
              //                                   SP.state[c_x].j_2, triP.J_2);
              if (part6_t == 0)
                continue;
              part6_t =
                  part6_t * PC.NinJ_quick(2 * jaco.l12, 2 * jaco.l3, 2 * L123,
                                          S12_2, 1, S_2, jaco.J12_2, jaco.j3_2,
                                          jaco.J_2);
              // part6_t = part6_t * ninJ.quick(2 * jaco.l12, 2 * jaco.l3,
              //                                2 * L123, S12_2, 1, S_2,
              //                                jaco.J12_2, jaco.j3_2,
              //                                jaco.J_2);
              // part6_t = part6_t * NineJSymbol_2(2 * jaco.l12, 2 * jaco.l3,
              //                                2 * L123, S12_2, 1, S_2,
              //                                jaco.J12_2, jaco.j3_2,
              //                                jaco.J_2);
              if (part6_t == 0)
                continue;
              phase_t = (S_2 + triP.J_2) / 2;
              part6_t = part6_t * (S_2 + 1) * phase(phase_t);

              part6 += part6_t;
            }
            part5 += part5_t * part6;
          } // nvtxRangePop();
          part4 += part4_t * part5;
        } // nvtxRangePop();
        part3 += part3_t * part4;
      } // nvtxRangePop();
      part2 += part2_t * part3;
    } // nvtxRangePop();
  }   // nvtxRangePop();
  result = part1 * part2;
  return result;
}

void TransCoeff::build_JTE(int J2, int T2, int E) {
  int ch_jac = jaco_nl.FindC(J2, T2);
  int ch_tri = TP.FindC(J2, T2);
  if (E > jaco_nl.channel[ch_jac].E_max || E < jaco_nl.channel[ch_jac].E_min) {
    cout << "wrong at TransCoeff::build_JTE" << endl;
    exit(0);
  }
  if (E > TP.channel[ch_tri].E_max || E < TP.channel[ch_tri].E_min) {
    cout << "wrong at TransCoeff::build_JTE" << endl;
    exit(0);
  }
  int beg_jac = 0;
  int beg_tri = 0;
  if (E > 0) {
    beg_jac = jaco_nl.channel[ch_jac].indexE[E - 1] + 1;
    beg_tri = TP.channel[ch_tri].indexE[E - 1] + 1;
    for (int i = E - 1; i > 0; i--) {
      if (jaco_nl.channel[ch_jac].indexE[i] > 0) {
        beg_jac = jaco_nl.channel[ch_jac].indexE[i] + 1;
        break;
      }
    }
    for (int i = E - 1; i > 0; i--) {
      if (TP.channel[ch_tri].indexE[i] > 0) {
        beg_tri = TP.channel[ch_tri].indexE[i] + 1;
        break;
      }
    }
  }
  int end_jac = jaco_nl.channel[ch_jac].indexE[E];
  int end_tri = TP.channel[ch_tri].indexE[E];
  // int end_jac = beg_jac + 10;
  // int end_tri = beg_tri + 10;
  if (-1 == 1) {
    cout << "J2 : " << J2 << " T2 : " << T2 << " E : " << E << endl;
    cout << "Jacobi_NL channel: " << ch_jac << "\t begin: " << beg_jac
         << "\t end: " << end_jac << endl;
    cout << "TP channel: " << ch_tri << "\t begin: " << beg_tri
         << "\t end: " << end_tri << endl;
  }
  int size_jac = end_jac - beg_jac + 1;
  int size_tri = end_tri - beg_tri + 1;
  Tval_JTE.resize(size_jac);
  for (int i = 0; i < size_jac; i++) {
    int i_x = i + beg_jac;
    Tval_JTE[i].resize(size_tri);
    for (int j = 0; j < size_tri; j++) {
      int j_x = j + beg_tri;
      Tval_JTE[i][j] =
          cal5(jaco_nl.state_NL[ch_jac][i_x], TP.state[ch_tri][j_x]);
    }
  }
  for (int i = 0; i < size_jac; i++) {
    int i_x = i + beg_jac;
    for (int j = 0; j < size_tri; j++) {
      int j_x = j + beg_tri;
      if (std::abs(Tval_JTE[i][j]) > 0.0000001) {
        cout << "\t" << i_x << "\t" << j_x << "\t" << Tval_JTE[i][j] << endl;
      }
    }
  }
}

void TransCoeff::build() {
  // Build Tval_vec[channel_abc][index_jaco_nl][index_abc]
  int cha_tri_num = TP.channel.size();
  Tval_vec.resize(cha_tri_num);
  Tval_vec_flag.resize(cha_tri_num);

  cout << "Begin to initial T-Coefficients " << endl;

  // int plan_t = 0;
  // int pos = 0;
  // int num = cha_tri_num - 1;
  Pgs_Bar pgs;
  pgs.init(cha_tri_num, 30);
  // nvtxRangePush("Trans_build : cha_tri loop");
  for (int cha_tri = 0; cha_tri < cha_tri_num; cha_tri++) {
    pgs.PrintBar(cha_tri);

    int J2 = TP.channel[cha_tri].J_2;
    int T2 = TP.channel[cha_tri].T_2;
    int cha_jac = jaco_nl.FindC(J2, T2);
    int Emin = TP.channel[cha_tri].E_min;
    int Emax = TP.channel[cha_tri].E_max;
    Tval_vec[cha_tri].resize(Emax + 1);
    Tval_vec_flag[cha_tri].resize(Emax + 1);

    // cout<<"Tri channel : "<<cha_tri<<"\t JacoNL channel : "<<cha_jac<<endl;
    // nvtxRangePush("Trans_build : E loop");
    for (int E = Emin; E <= Emax; E++) {
      if (E > jaco_nl.channel[cha_jac].E_max ||
          E < jaco_nl.channel[cha_jac].E_min) {

        //  cout<<"wrong at TransCoeff::build(); Jaco_NL space shuold become
        //  bigger"<<endl;
        //  exit(0);
      }
      int beg_jac = jaco_nl.channel[cha_jac].indexE_beg[E];
      int end_jac = jaco_nl.channel[cha_jac].indexE[E];

      int beg_tri = TP.channel[cha_tri].indexE_beg[E];
      int end_tri = TP.channel[cha_tri].indexE[E];

      int tri_size = end_tri - beg_tri + 1;
      int jac_size = end_jac - beg_jac + 1;

      if (end_tri == 0 && beg_tri == 0 && TP.state[cha_tri][end_tri].E != E)
        tri_size = 0;

      Tval_vec[cha_tri][E].resize(jac_size);
      Tval_vec_flag[cha_tri][E].resize(jac_size);

      // cout<<" E : "<<E<<"\t beg_jac : "<<beg_jac<<" end_jac : "<<end_jac<<"\t
      // beg_tri : "<<beg_tri;
      // cout<<" end_tri : "<<end_tri<<endl;
      // nvtxRangePush("Trans_build : omp parrlel jac loop ");
      //#pragma omp parallel for schedule(dynamic, 1)
      for (int i = beg_jac; i <= end_jac; i++) {
        int i_index = i - beg_jac;
        Tval_vec[cha_tri][E][i_index].resize(tri_size);
        Tval_vec_flag[cha_tri][E][i_index].resize(tri_size);
        if (tri_size == 0)
          continue;

        for (int j = beg_tri; j <= end_tri; j++) {
          int j_index = j - beg_tri;
          Tval_vec[cha_tri][E][i_index][j_index] = 0.0;
          Tval_vec_flag[cha_tri][E][i_index][j_index] = 0;
        }
        // State_Jacobi Jaco_no_t = jaco_nl.state_NL[cha_jac][i];
        // nvtxRangePush("Trans_build :  tri loop ");
        //#pragma acc parallel loop copyin(TP.state[cha_tri][beg_tri:end_tri],
        // Jaco_no_t) copy(Tval_vec_t, Tval_vec_flag_t) #pragma acc parallel
        // loop
        for (int j = beg_tri; j <= end_tri; j++) {
          if (TP.state[cha_tri][j].F < 0 && TP.state[cha_tri][j].I < 0)
            continue;
          int j_index = j - beg_tri;
          // double result =
          //     cal5_OA(Jaco_no_t, TP.state[cha_tri][j]);
          double result =
              cal5(jaco_nl.state_NL[cha_jac][i], TP.state[cha_tri][j]);
          //
          Tval_vec[cha_tri][E][i_index][j_index] = result;
          Tval_vec_flag[cha_tri][E][i_index][j_index] = 1;
        }
        // nvtxRangePop();
      }
      // nvtxRangePop();
    }
    // nvtxRangePop();
  }
  // nvtxRangePop();
  pgs.EndBar();
  cout << " T-Coefficient initial finished " << endl;
}

void TransCoeff::build_OA() {
  // Build Tval_vec[channel_abc][index_jaco_nl][index_abc]
  int cha_tri_num = TP.channel.size();
  Tval_vec.resize(cha_tri_num);
  Tval_vec_flag.resize(cha_tri_num);

  cout << "Begin to initial T-Coefficients " << endl;
  Binom binom_x;
  int n_max_x = 100;
  binom_x.init(n_max_x);

  CG_coeff cg_x;
  NinJ_coeff ninJ_x;
  //
  // Hobra hobra_xt(binom_x, cg_x, ninJ_x);
  // // hobra_x
  // hobra_xt.init_2(binom_x, cg_x, ninJ_x);
  // int plan_t = 0;
  // int pos = 0;
  // int num = cha_tri_num - 1;

  CG_coeff cg_t;
  cg_x.init(15);

  SixJ_coeff sixJ_t;
  sixJ_t.init(14);

  Pre_coeff PC;
  PC.init_Binom();
  PC.get_CG(cg_t);
  PC.get_SixJ(sixJ_t);
  cout << " Finished Init PC " << endl;
  Time_Count TC_t;
  TC_t.start();
  // int size_vec_sixJ = 0;
  // double binom_sum_t = 0;
  // for (int i = 0; i < 30; i++) {
  //   binom_sum_t += hobra_x.fac_cal(i);
  // }
  // double binom_sum_t_OA = 0;
  // #pragma acc parallel loop reduction(+ : binom_sum_t)
  //   for (int i = 0; i < 30; i++) {
  //     binom_sum_t_OA += PC.fac_cal(i);
  //     size_vec_sixJ = PC.vec_6j.size();
  //   }
  //   cout << " !!! binom_sum_t = " << binom_sum_t
  //        << "\t binom_sum_t_OA = " << binom_sum_t_OA << endl;
  //   cout << "\t SixJ vector size : " << size_vec_sixJ << endl;
  // exit(0);

  int size_sp = SP.state.size();
  std::vector<State_SP> SP_state(size_sp);
  for (int k = 0; k < size_sp; k++) {
    SP_state[k] = SP.state[k];
  }

  Pgs_Bar pgs;
  pgs.init(cha_tri_num, 30);

  nvtxRangePush("data copyin(SP_state)");
#pragma acc data copyin(SP_state, PC)
  {
    nvtxRangePop();
    nvtxRangePush("Trans_build : cha_tri loop");
    for (int cha_tri = 0; cha_tri < cha_tri_num; cha_tri++) {
      pgs.PrintBar(cha_tri);

      int J2 = TP.channel[cha_tri].J_2;
      int T2 = TP.channel[cha_tri].T_2;
      int cha_jac = jaco_nl.FindC(J2, T2);
      int Emin = TP.channel[cha_tri].E_min;
      int Emax = TP.channel[cha_tri].E_max;
      Tval_vec[cha_tri].resize(Emax + 1);
      Tval_vec_flag[cha_tri].resize(Emax + 1);

      vector<State_TP> TP_ch_t(cha_tri_num);
      int size_temp = TP.state[cha_tri].size();
      TP_ch_t.resize(size_temp);
      for (int i = 0; i < size_temp; i++) {
        TP_ch_t[i] = TP.state[cha_tri][i];
      }

      // cout<<"Tri channel : "<<cha_tri<<"\t JacoNL channel : "<<cha_jac<<endl;
      nvtxRangePush("Trans_build : E loop");
      for (int E = Emin; E <= Emax; E++) {
        if (E > jaco_nl.channel[cha_jac].E_max ||
            E < jaco_nl.channel[cha_jac].E_min) {

          //  cout<<"wrong at TransCoeff::build(); Jaco_NL space shuold become
          //  bigger"<<endl;
          //  exit(0);
        }
        int beg_jac = jaco_nl.channel[cha_jac].indexE_beg[E];
        int end_jac = jaco_nl.channel[cha_jac].indexE[E];

        int beg_tri = TP.channel[cha_tri].indexE_beg[E];
        int end_tri = TP.channel[cha_tri].indexE[E];

        int tri_size = end_tri - beg_tri + 1;
        int jac_size = end_jac - beg_jac + 1;

        if (end_tri == 0 && beg_tri == 0 && TP.state[cha_tri][end_tri].E != E)
          tri_size = 0;

        Tval_vec[cha_tri][E].resize(jac_size);
        Tval_vec_flag[cha_tri][E].resize(jac_size);

        // cout<<" E : "<<E<<"\t beg_jac : "<<beg_jac<<" end_jac :
        // "<<end_jac<<"\t beg_tri : "<<beg_tri; cout<<" end_tri :
        // "<<end_tri<<endl;
        nvtxRangePush("Trans_build : omp parrlel jac loop ");
        //#pragma omp parallel for schedule(dynamic, 1)
        vector<vector<double>> Tval_vec_t(jac_size);
        vector<vector<int>> Tval_vec_flag_t(jac_size);

        for (int i = beg_jac; i <= end_jac; i++) {
          int i_index = i - beg_jac;
          Tval_vec[cha_tri][E][i_index].resize(tri_size);
          Tval_vec_flag[cha_tri][E][i_index].resize(tri_size);
          Tval_vec_t[i_index].resize(tri_size);
          Tval_vec_flag_t[i_index].resize(tri_size);
          if (tri_size == 0)
            continue;

          for (int j = beg_tri; j <= end_tri; j++) {
            int j_index = j - beg_tri;
            Tval_vec[cha_tri][E][i_index][j_index] = 0.0;
            Tval_vec_flag[cha_tri][E][i_index][j_index] = 0;
            Tval_vec_t[i_index][j_index] = 0.0;
            Tval_vec_flag_t[i_index][j_index] = 0;
          }
        }

        vector<State_Jacobi> Jaco_no_t(jac_size);
        for (int i = 0; i < jac_size; i++) {
          Jaco_no_t = jaco_nl.state_NL[cha_jac];
        }
        nvtxRangePush(" loop for cal5_OA ");
#pragma acc data present(SP_state)
#pragma acc parallel copyin(Jaco_no_t, TP_ch_t) copy(Tval_vec_t, Tval_vec_flag_t) vector_length(32)
        {
#pragma acc loop gang independent collapse(2)
          for (int i = beg_jac; i <= end_jac; i++) {
            for (int j = beg_tri; j <= end_tri; j++) {
              // if (TP.state[cha_tri][j].F < 0 && TP.state[cha_tri][j].I < 0)
              //   continue;
              int i_index = i - beg_jac;
              if (tri_size == 0)
                continue;
              if (TP_ch_t[j].F < 0 && TP_ch_t[j].I < 0)
                continue;
              int j_index = j - beg_tri;
              double result = cal5_OA(Jaco_no_t[i], TP_ch_t[j], SP_state, PC);
              Tval_vec_t[i_index][j_index] = result;
              Tval_vec_flag_t[i_index][j_index] = 1;
            }
          }
        }
        nvtxRangePop();
        nvtxRangePush(" Grep Tval_vec[cha_tri][E][i_index][j_index] ");
        for (int i = beg_jac; i <= end_jac; i++) {
          int i_index = i - beg_jac;
          for (int j = beg_tri; j <= end_tri; j++) {
            int j_index = j - beg_tri;
            Tval_vec[cha_tri][E][i_index][j_index] =
                Tval_vec_t[i_index][j_index];
            Tval_vec_flag[cha_tri][E][i_index][j_index] =
                Tval_vec_flag_t[i_index][j_index];
          }
        }
      } // E
      nvtxRangePop();
    } // ch_tri
    nvtxRangePop();
  }

  nvtxRangePop();
  pgs.EndBar();
  cout << " T-Coefficient initial finished " << endl;
  cout << " T-coefficients Time cost without Init PC: " << TC_t.sec() << "s"
       << endl;
}

void TransCoeff::print(std::string filename) {
  std::ofstream file;
  file.open(filename);
  int size_ch = Tval_vec.size();
  for (int ch = 0; ch < size_ch; ch++) {
    file << "Channel TP : " << ch << endl;
    int size_E = Tval_vec[ch].size();
    for (int E = 0; E < size_E; E++) {
      file << "=== E : " << E << "  ===" << endl;
      file << "JacoNL index"
           << "\t TP index" << endl;
      int size_jac = Tval_vec[ch][E].size();
      for (int i = 0; i < size_jac; i++) {
        int size_TP = Tval_vec[ch][E][i].size();
        for (int j = 0; j < size_TP; j++) {
          if (std::abs(Tval_vec[ch][E][i][j]) > 0.00000001)
            file << i << "\t\t" << j << "\t\t" << Tval_vec[ch][E][i][j] << endl;
        }
      }
    }
  }
  file.close();
}
